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Abstract 

We compare the DAMA signal to the muon flux underground, which is annually mod- 
ulated due to temperature variations in the stratosphere. We show that the muon flux at 
LNGS and the DAMA signal are tightly correlated. Different mechanisms were proposed in 
the literature by which muon-induced events may dominate the signal region in DAMA. We 
discuss simple statistical constraints on such mechanisms and show that the DAMA collab- 
oration can falsify the muon hypothesis, if it is wrong, by reporting their annual baseline 
count rates. 

1 Introduction 

The DAMA/LIBRA experiment [1] is a radiopure Nal(Tl) scintillation detector located at the 
Gran Sasso National Laboratory (LNGS), searching for the annual modulation signature of dark 
matter particles. With integrated exposure larger than 1 ton x year and covering 13 annual cycles, 
DAMA/LIBRA and the former DAMA/Nal results reveal a clear annual modulation in the count 
rate of scintillation events [21 [3]. The modulation period, phase and amplitude and the fact that 
modulation is only reported to exist in single-hit events, where only one Nal(Tl) crystal out of 
the 25-crystal array is triggered, are all consistent with expectations from direct detection of 
dark matter [1]. 

A thorough investigation of possible background processes which may explain the DAMA 
result is clearly of utmost priority [SJ [HI [Zl E] ■ The DAMA collaboration argues that no back- 
ground was found which satisfies all the features attributed to the modulation signal |2j , including 
the time dependence, amplitude and event distribution in the detector array. Despite of those 
arguments, one particular source of background, common to all underground low-noise experi- 
ments, merits very careful consideration when one tries to search for an annual modulation. This 
background is the flux of penetrating underground muons [5J [7] . 

Annual variations in the muon flux underground have been measured by many experiments 
ITOl [TTl [T2l fT3] . Modulation was also reported for neutrons [14] which are produced, among other 
mechanisms, through muon interactions with rock. In the northern hemisphere the muon rate 
peaks around June- July, close to the DAMA phase. The coincidence of the modulation in phase 
and amplitude motivated Refs. [BJ [7^ to consider underground muons and/or neutrons as a 
possible explanation to the DAMA signal. 

In this paper we investigate further the hypothesis that muons may explain the DAMA 
anomaly and propose simple methods to test it. In Sec. [2] we extend the discussion in [6j [7] 
by comparing the DAMA signal to the muon background directly. We point out that a muon 
explanation for DAMA cannot rely on muons directly traversing the detector material, but 
instead should involve muon interactions in a larger volume surrounding the detector. We study 
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the time behavior of the muon rate and show that the DAMA coUaboration can rule out the 
muon hypothesis, if it is wrong, using only existing data. This can be done exploiting the long 
term variability of the baseline muon flux, which is not expected from a dark matter signal. We 
conclude in Sec.[3l In App.[2]we recapitulate the physics of the muon time variability. 



2 Testing the muon background hypothesis 

We begin by describing the DAMA signal. DAMA applies a software threshold of 2 KeV, in 
order to avoid the tail of PMT noise peaking at lower energies. Here, for concreteness, we focus 
on the energy interval 2-4 KeV. In this interval, using Refs. [151 [5] we obtain the baseline count 
rates averaged over four DAMA/LIBRA and two DAMA/Nal annual cycles, respectively: 



s = 1.15 cpd/kg/KeV (DAMA/LIBRA 1-4), 
s = 1.45 cpd/kg/KeV (DAMA/Nal 3 - 4). 



(1) 



While Eq. ([T]) gives the baseline rate averaged over a few cycles, it is important to note that 
DAMA defines the "baseline count rate" separately for each annual cycle 0. The residual count 
rate Ss in DAMA/Nal and DAMA/LIBRA was given in [2] and [5]. The residual count rate in 
each annual cycle is obtained by subtracting the baseline rate of that cycle from the data. 

Next, we require information on the muon flux. Measurements of the daily muon intensity 
at LNGS were carried out by the LVD experiment [5]. In Fig. [T] we plot the DAMA/LIBRA and 
DAMA/Nal residuals (blue), in the form of percent modulation relative to the mean baseline 
count rates of Eq. ^ , on top of the LVD muon intensity modulation (green) In presenting 
the muon modulation we measure the residuals with respect to the baseline muon rates at 
each annual cycle separately, starting each cycle on Sep 9. Fig. [1] makes clear that the muon 
background should be regarded with care. 



(0 

■O -2 



DAMA Nal 



I 



,1 1 






t 

















ICARUS n 




2500 3000 3500 

days since 1 Jan 1995 



Figure 1: DAMA residuals (blue [H [3]) and LVD muon intensity residuals (green [3]) in per- 
cent from the respective baselines. ICARUS neutron measurements during 1997-1998 are added 
(red [H]). 



^We thank Rita Bernabei for clarification on this point. 

^The pubhcly available LVD data [9] is provided in a dense graphical form. In converting the muon data to 
digital form, we end up with a reduced data set sampled roughly on a five-day basis. Our procedure carries the 
LVD daily error bars to the digitized five-day data, losing a potential factor of two improvement in accuracy. 
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The DAMA collaboration considered the possibility that their signal is associated with 
muons [11] or with the fast neutron flux that muons induce [TTj H]. These possibilities were 
discarded based on three types of arguments: 

1. Intensity: Considering fast neutrons, a rough estimate of the flux was made, indicating 
that the fast neutron flux is of order 10^^ cpd/kg/KeV, three orders of magnitude smaller 
than the DAMA count rate. This small flux was discarded as harmless. 

2. Background rejection: Muon-induced processes (including fast neutrons) were argued to 
fail a number of background rejection criteria. The relevant criteria include the presence 
of multiple-hit events and deposition of energy above the signal region. 

3. Modulation phase: The authors of [TB] entertained the possibility that muon-induced events 
evade somehow the background rejection and intensity arguments, and considered the muon 
flux modulation of LVD [3]. Based on sinusoidal fits to both observables - the DAMA 
residuals and the muon rate - Ref. [16] concluded that the inferred phases of the two 
effects are separated by more than Str. 

All of the three arguments above were challenged in ^ and (71. 

First, two mechanisms were discussed wherein muon-induced interactions in and around the 
DAMA apparatus could evade both the intensity and the background rejection criteria. Ref. [5] 
discusses in detail the neutron background, arguing that it is possible for neutrons to deposit 
energy in the DAMA signal region without impacting the higher energy region, and showing that 
the yield of muon-induced neutrons may have been underestimated by orders of magnitude. The 
DAMA signal could be associated to ~ 3 KeV Auger-L electrons and X-rays accompanying the 
decay of (lifetime ~ 36 minutes), formed from ^^''l by neutron capture. Ref. [7] points out to 
the possibility of delayed phosphorescence in Nal(Tl), which could be initiated by muon-induced 
energy deposition inside the detector but trigger long after the actual muon passage, evading 
coincidence veto. We have nothing to add to the discussions of background rejection in [6] [7], 
and we refer the reader to these references for more details. We do have further insight regarding 
the intensity argument, which we consider in Sec. 12.11 

Second, the phase argument was challenged. Here, too, our analysis adds significant infor- 
mation. We study the phase argument in more detail in Sec. 12.21 

2.1 Direct and indirect muons 

Without referring to the precise nature of muon-induced events (in particular, whether they 
involve secondary neutrons or not), there are two basic ways by which muons could induce the 
DAMA signal: directly, following muon passage through the detector, or indirectly, e.g. by muon 
interactions in a larger volume surrounding the detector. We now show that current LVD and 
DAMA data suffice to discriminate between these two possibilities, making the first one highly 
unlikely, regardless of the detailed physics responsible for converting muon energy deposit to 
single hit readouts. 

First, one may suppose that DAMA counts arise from muons traversing the Nal(Tl) detector 
material itself [7], or producing an enhanced yield of neutrons in the 15cm-thick led shield 
attached to the detector [B]. The rate of direct muons in the DAM A/LIBRA detector is 

-Rdiroct^ = A X w 10 muons/day, (2) 

where we used the LVD intensity measurement « 28 muons/m^/day and estimated the effec- 
tive area of the DAMA/LIBRA detector by A « 0.35 m^. For comparison, the DAMA/LIBRA 
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count rate in the 2-4 KeV signal region is s M AE sa 530 cpd, as can be read from Eq. ([T]) using 
M =232 kg for the active scintillator mass relevant for the first five annual cycles [2- 

We learn that the direct muon possibility requires a large yield, y « 50 signal counts/muon, 
to accommodate the signal rate. By itself, the requirement of a large yield does not seem 
particularly constraining. As noted in [7], muons traversing the DAMA apparatus deposit some 
©(GeV) energy in the detector, leaving room for small efficiencies of O (10~*) in producing tens 
of KeV light pulses. The problem with direct muons is not related to efficiencies. Rather, it has 
to do with the fact that a small number of seed muons are taken to be responsible for a much 
larger number of eventual counts. 

Regardless of the statistical distribution of the yield y, the relative spread of the induced 
DAMA counts in a time bin ti cannot be smaller than that of the small number of direct muon 
events: 



> 



-^direct fj,,i 



30% X t-'/\ (3) 



with ti given in days, Si the DAMA count rate per unit mass per unit energy and Cs. the error 
bar on the rate (assumed to be statistics dominated). For a typical integration time of order 50 
days, Eq. ([3]) predicts a spread of at least four percent in DAMA/LIBRA residuals, a factor > 5 
larger than the spread visible in practice. We conclude that the direct muon scenario is unlikely. 
This scenario would predict a much larger spread to the DAMA residuals than that observed Q 
The second possibility is that DAMA events are induced by muon interactions in a large 
volume around the detector [6l [HI [19]. In this case, the yield y can be 0(1). For given yield y 
we can estimate the effective area through which muons should be collected, 

DAMA/LIBRA counts/day s M AE 20 , 

Acff = / 2,A = f ~ — "1 • 4 

y X muons/m^/day yl^ y 

In Eq. we took muons to cause the signal in the 2-4 KeV energy range with AE = 2 KeV, 
where DAMA observe most of the modulation. If muons are responsible to the signal in a larger 
energy range, then A^g should grow proportionally. 

Note that the mean single hit rate reported by DAMA is essentially flat as a function of energy 
at least up to 10 KeV, except for a shallow bump of order 10 percent of the baseline which occurs 
between 2-4 KeV0. As pointed out in [7], demanding the muon signal to diminish below 10 KeV 
implies that some other (unmodulated) background must grow to keep the observed spectrum 
flat. As also pointed out in [7], this curious fine tuning problem pertains to the dark matter 
hypothesis just as well as it does to background-related alternatives. 

We can use the DAMA data to constrain the yield y. To do this, we compute the quantity 



M,AEe,U, (5) 



where Mi is the exposed mass and is the duty cycle during the time bin ti. If the DAMA 
signal originates from a seed Poisson process with average yield y, then one expects (qi) « y. 
Using the 2-4 KeV DAMA/LIBRAl-5 residuals data with Mi = 232 kg and = 60%, we obtain 
(q) ^ 1.7. After comparing the distribution of qi obtained from the true data to that from 



^Eq. ^ assumes that all of the y ~ 50 counts arising from a direct muon hit occur in the time bin ti containing 
the hit. This reasoning is modified if the counts exhibit very long delay tdolay ^ after the muon passage, in 
which case we would need to replace ti — > tdclay in Eq. (|3jl. However, to remedy the discrepancy in direct muon 
statistics, tdelay ~ 20 > 1 year is required, which (a) seems unrealistically long and (b) implies a very extended 
rise time to the count rate, absent in the data. 

■^The location of the 2-4 KeV bump agrees with expectations from internal '^'^K contamination \5\. 
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several sets of cosine-modulated Poisson random processes, we conclude that any background- 
related explanation for the DAM A signal should exhibit y < 3. 

Using Eq. with y — 1.7, in agreement with the scaling of DAMA/LIBRA error bars, 
gives Acs ~ 12 m^. This could imply that the relevant muons interact within a volume extend- 
ing roughly 1.5-2 meters to either side from the center of the detector, the dimensions of the 
DAMA/LIBRA apparatus including the surrounding meter-thick concrete housing [I]. 

2.2 Temporal correlation 

Ref. [5] compared the DAMA signal with ICARUS measurements of the neutron flux at LNGS [U] , 
finding a troubling correlation. We show the ICARUS data in Fig. [T] (red marks). The ICARUS 
data set contains five measurements, extending from April 1997 to May 1998, in which neutron 
events are associated with proton recoils above 3.5 MeV in a home-made organic liquid scintil- 
lation detector. Clear modulation of order few percent is present. Ref. [B] fitted a sinusoidal 
curve to the ICARUS data, discovering that the same fit interpolates nicely the DAMA residuals. 
While certainly alarming, the correlation reported in [6] is based on extrapolation of a loosely 
constrained fit. 

Ref. [7| argued against the statistical inference employed by the DAMA collaboration in [IB] 
to reject the muon hypothesis, making the following straightforward observations: 

• In quoting a > 5(t discrepancy, DAMA ignored the reported error on the LVD muon fit. 
In fact, the reported LVD phase is July 5±15 days. The combined DAMA/LIBRA and 
DAMA/Nal phase. May 26±7, is thus less than 3a away from the muon phase. 

• DAMA assumed as prior a one year period to their fit. Adding the period to the fit would 
further reduce the statistical discrepancy. 

• Importantly, Ref. |7j observe that the fine-sampled LVD muon data does not follow a true 
sinusoidal behavior. If the underlying process for DAMA is also non-sinusoidal, fits to the 
DAMA data in which a sinusoidal form was assumed as prior will be affected by systematics, 
potentially invalidating statistical inference. 

To explore the muon hypothesis in more detail, in Fig. [2] we bin the muon residuals on the 
same binning as the DAMA residuals. To emphasize the time scale, a 40 day period - equal to 
the difference between LVD's and DAMA's sinusoidal best fit phases - is depicted in Fig. [5] as a 
black line above day 3500. 

To estimate the significance of a phase discrepancy between the DAMA and muon data we 
proceed as follows. We define the muon hypothesis as saying that the DAMA count rate in time 
bin ti is given by 

where N^^i is Poisson distributed with mean 

{Nf_,,i) = Acs If_,,ieiti, (7) 

with Acff taken from Eq. Q and = 60%. In Eq. ([7]), /^.i is the muon intensity in time 
bin ti . For simplicity, we adopt the yield y = 2 and neglect sources of background other than 
those correlated with muons. We generate multiple realizations of N^^i for the time bins of 
DAMA/LIBRAl-5. For each realization we define and remove the annual baselines, and then 
fit the resulting time sequence to a cosine function, mimicking the DAMA procedure. We thus 
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Figure 2: DAM A residuals (blue) and binned muon intensity residuals (green). ICARUS neutron 
measurements during 1997-1998 are added (red). 



obtain distributions of the following quantities: (1) the best fit phase to, allowing both the 
amplitude and the period to float; (2) the best fit x^, minimizing over the amplitude with period 
and phase fixed to 1 year and day 152.5, respectively. 

The resulting distributions of ^-^d are shown in the left and right panels of Fig. [3] 
(green). The number of dof in the fit are 37, with 38 DAMA/LIBRAl-5 time bins and one 
fitted amplitude parameter. What interests us here is the distribution of and to for the mock 
data, in comparison with the results x^ jdo] « 1 and to ~ 140 days [3], obtained by DAM A 
for their true data and shown as blue bands in Fig. [31 We conclude that the simplistic muon 
hypothesis defined by Eqs. (|6][7l) is consistent with the results reported by DAMA. 




Figure 3: Distributions of the best fit (left panel) and phase (right panel), generated for 
multiple realizations of mock data under the muon hypothesis. DAMA fit results are shown as 
blue vertical bars. 
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2.3 The usefulness of annual baseline data 



If DAMA counts arise from muons - be it through secondary particles produced by interactions 
in the surrounding volume, activation of impurities or even in some mundane manner via distur- 
bances to the electronics - then the yearly baseline count rates should reflect long term trends 
in the muon data. In the energy range 2-4 KeV, these annual baseline rates should be measured 
to accuracy of ~ 0.25%. 

In Fig. 0] we plot again the DAMA and muon data. This time, instead of removing the 
baselines for each year separately, we show the muon residuals defined with respect to the seven 
year average in the period 9 Sep 2001 until 9 Sep 2008. The muon intensity exhibits long term 
modulation with period ~ 6 years and non negligible amplitude > 1%. The annual baselines 
of the muon intensity are shown as dark green bands. The width of each band is ±0.25%, the 
estimated DAMA/LIBRA resolution for a one year integration time. 




500 1000 1500 2000 2500 3000 3500 4000 4500 5000 



days since 1 Jan 1995 

Figure 4: Long term modulation in the annually-averaged muon rate (thick dark green) with 
respect to the seven year average in the period Sep 9 2001 until Sep 9 2008. The thickness of the 
horizontal bars is ±0.25%, representing the expected DAMA resolution for the energy range 2-4 
KeV. Here, in contrast to Fig. [TJ the short term muon intensity residuals (green) were defined 
with respect to the seven year baseline rather than to the different annual baselines. 



3 Conclusions 

The DAMA anomaly has been in existence for over a decade. No solution in terms of a back- 
ground process has thus far been verified. In this paper we compared the DAMA signal to the 
underground muon flux at LNGS. We found that the two observables are tightly correlated. 
The correlation persists throughout the 13 annual cycles of operation of DAMA/LIBRA and its 
former DAMA/Nal (see Figs. [1] and |6]), motivating muon-induced events as an explanation for 
the DAMA result. 

Our work adds to Refs. which suggested various mechanisms by which muon-induced 

interactions could potentially bypass the DAMA background rejection criteria, and pointed out 
that both the amplitude and the phase of the annual modulation in the underground flux of 
neutrons and muons at LNGS are not inconsistent with DAMA's. 

The underground muon flux exhibits short and long term variations of order percent on 
time scales ranging from weeks to several years, deviating systematically from pure sinusoidal 
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form. This behavior is not expected from the dark matter signal The DAMA cohaboration can 
potentiahy exclude a muon-induced background by reporting their high statistics annual baseline 
count rates in the signal region. 

While there are suggestions, the actual mechanism by which muons may eventually lead to 
the DAMA signal remains unknown. Any hints which may hide in the data are thus important. 
Inspecting the spread of DAMA residuals, we showed that it is unlikely for the DAMA signal to 
arise solely from muons traversing the detector material itself. This scenario, suggested in [7], 
would lead to a large irreducible statistical spread in DAMA events, which is not manifest in the 
data. In general, the same argument disfavors high yield processes in which few seed interactions 
deliver a large amount of energy into the DAMA detector material, subsequently transformed 
into multiple KeV light pulses in the signal region. Using available DAMA data we showed that 
the yield of such process must be moderate, y <3. The signal may arise from muon interactions 
throughout a volume extending 2 — 3 meters to either side of the active detector. 

Lastly, when considering the implications of Fig. [TJ a comment is in order. One should bear 
in mind the accidental phase proximity of the expected dark matter peak (early June) and solar 
solstice (late June, around which time the atmosphere in the northern hemisphere is warmest 
leading to a maximum in the underground muon intensity). This accidental coincidence should 
be regarded as a prior, making it harder to discriminate a muon-induced background from the 
dark matter hypothesis [5D], but not adding to the relative significance of the former. Having 
made this cautionary comment, it is still remarkable that the modulation amplitudes of the 
DAMA signal and the muon intensity are so close. 



A The muon flux — temperature correlation 

High energy muons are generated through meson decays in the stratosphere and penetrate deep 
underground. Annual temperature changes induce a corresponding modulation in the muon rate, 

f- = o.^. (8) 

Hp. leff 

where ax is a proportionality coefficient which depends on the cosmic ray composition in the 
upper levels of the atmosphere, and T^ff is an effective temperature [211 llOj . Here we provide a 
quick analytical derivation of Eq. ([5]) using a toy model. We note that Eq. (jS]) is well known in the 
literature since early times [221 [531 HU [TU] ^nd so this part is mainly added for completeness. We 
then proceed to compute Tg/ / and use it first as independent verification of the time variability 
of the muon rate on different time scales and, second, to cover the entire 13 years of DAMA/Nal 
and DAMA/LIBRA operation. 



Understanding the muon flux temperature correlation. Let us compute the effect of 
temperature variations on the muon intensity in a toy model, approximating the atmosphere as 
a single layer of air with uniform density and temperature, considering an isotropic distribution 
of muons, and ignoring the precise cosmic ray composition and the details of muon propagation 
in rock. We work with energy in GeV units, e = E/GeV. 

Atmospheric muons are produced in the decay of pions and kaons, generated in turn by 
cosmic ray collisions with air. Here we consider only the dominant pion source. At production, 
pions inherit the spectral index of the primary cosmic rays such that their flux obeys 

J.(e) K (9) 
with 7 = 2.75. The spectrum of pion-decay muons depends on the competition of two processes: 
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1. pion decay, controlled by the observer frame lifetime t7r(e) = -^p^ = 1.8 -10 e s 

2. collisions with air, with characteristic time tc — ^ 10^^ ^ s, the air density n given 
in units of 3 • 10^^ cm~^, the relevant order of magnitude at the tens of km height where 
most of the muons form. The air density responds to temperature changes, — o^T . Thus 
tc oc T. 

Equating the time scales i^r and tc defines a single critical energy: 

t^{t) = tc^e='^^im-. (10) 

At energy e collision losses and decay are equally important. We see that e oc T. 

Collisions with air degrade the energy of pions, such that the differential probability of a pion 
produced at energy e-j^^Q to decay at energy e-^^d is: 



We also need the distribution of muons from pion decay: 



(11) 



[ , o.w. J 

Using Eqs. (|lim2p we have 

■Ju,\^ii) — I <l^TT,dJii{^-ir.d,^n) I «£;!•, 0J7r(,e7r,0j ; — "j \ ~ 



X Jr^ Jo 



7-2g-irKl-?)_ (13) 



For X < 0.1, T{x) « 0.3. For x > 1, T{x) « 0.13/2;. Thus for < e, the muon flux at sea level 
inherits the production spectrum of the pions, while for > e the muon flux is softer by one 
power of energy. 

Computing the muon rate above some threshold eth gives: 

/"OO /"OO 

i?oc / deJ^(e) cx e^-'^ / dxx"^I'(x). (14) 



We learn that if tth > e, then R oc e, while if eth < e, then R is approximately independent of e. 
As shown above, e oc T and thus for eth > e we find 

AR AT 

-^ = -T — , (15) 

with aT ~ 1. This completes the derivation of Eq. ([5]), up to refined considerations which are 
needed in order to determine the precise value of a-r- Note that a 3 kmwe rock overburden 
at LNGS is equivalent to eth > 10'^ (1 TeV) for the surface muon energy, sufficiently large 
compared to e for pions. The MACRO collaboration [T^ measured the value of ar at their 
apparatus in LNGS, obtaining = 0.98 ±0.12 for their full data set in the years 1991-1994 and 
ar = 0.83 ± 0.13 for a partial data set extending over 1993-1994. 
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Computation of Teff. The effective temperature Te// corresponds to an average of the tem- 
perature at different layers in the atmosphere, weighted by the probability for the observed muons 
to be formed at each layer. A useful approximate relation is given by |10j 

where Ajr = 160 gcm^^, A^r — 120 gcm^^ and the integration is over the slant depth X (in units 
of gcm^^), for which we use an approximate isothermal expression |21] 

X{h) ^ Xoe-''/''", Ao = 1030 gcm~^ /lo = 6.4 km (17) 

with h the geophysical height in km. 

We obtain atmospheric temperature records from the Integrated Global Radiosonde Archive 
(IGRA) of the National Climatic Data Center (NCDC), US Department of Commerce [M] . 
Archival information is available for the Pratica di Mare station in Italy, located at a distance 
of 124 km from LNGS. IGRA data contains temperature measurements at different discrete 
geophysical heights, typically derived up to few times per day. Using these data we compute 
Teff, approximating the integrals in Eq. (|T6)) by discrete sums. 

In Fig. [5] we plot the temperature modulation above the Pratica di Mare station (green) 
together with the LVD muon rate data (black markers with red sine fit, adopted from [S]). The 
correlation is clear. Occasional systematic deviations can be traced to the geographical distance 
between the weather station and LNGS. 




Figure 5: Muon flux at LVD (black data, red sinusoidal fit, taken from [9 ) and ar (AT/T) at 
the stratosphere (green). 



Lastly, in Fig. [B] we plot the DAMA residuals together with the daily stratospheric temper- 
ature modulation. The temperature analysis allows us to follow the 13 years of DAMA data 
taking, extending beyond the 8 years of LVD muon records. 
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days since 1 Jan 1995 



Figure 6: DAMA residuals (blue) and stratospheric temperature residuals /ST^f f /T^f f (green), 
in percent from the respective baselines. 
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